big_func3 <- function(s){

  if(s == 0){
    bs_s <- 1:n
  } else {
    bs_s <- sample(1:n,n,TRUE)
  }

  # Compute bootstrap sample equivalent of vectors
  ids <- id[bs_s]
  y0s <- y0[bs_s]
  y1s <- y1[bs_s]
  y2s <- y2[bs_s]
  y3s <- y3[bs_s]
  X1s <- X1[bs_s,]
  X2s <- X2[bs_s,]
  X3s <- X3[bs_s,]

  X1_wide_s <- X1_wide[bs_s,]
  X2_wide_s <- X2_wide[bs_s,]
  X3_wide_s <- X3_wide[bs_s,]

  # Additional definitions for particular models.
  y_lm <- c(y1s,y2s)
  X_wide_cs <- rbind(X1_wide_s,X2_wide_s)

  ## Linear
  lin1b <- lm(y_lm ~ X_wide_cs)
  # Identical to:
  # lm(health~child+married+emp, df_country %>% filter(period %in% c(2,3)))

  r_lin1b <- tibble(
    estimator = "LM-WIDE",
    s = s,
    point = ifelse(s==0,"point","bs_rep"),
    b_ltdhi = lin1b$coefficients[2],
    b_child = lin1b$coefficients[3],
    b_married = lin1b$coefficients[4],
    b_male = lin1b$coefficients[5],
    b_emp_2 = lin1b$coefficients[6],
    b_emp_3 = lin1b$coefficients[7],
    b_emp_4 = lin1b$coefficients[8],
    b_urbanisation_2 = lin1b$coefficients[9],
    b_urbanisation_3 = lin1b$coefficients[10],
    b_ageg_2 = lin1b$coefficients[11],
    b_ageg_3 = lin1b$coefficients[12],
    b_ageg_4 = lin1b$coefficients[13],
    b_ageg_5 = lin1b$coefficients[14],
    b_ageg_6 = lin1b$coefficients[15],
    b_educ6_1 = lin1b$coefficients[16],
    b_educ6_2 = lin1b$coefficients[17],
    b_educ6_3 = lin1b$coefficients[18],
    b_educ6_4 = lin1b$coefficients[19],
    b_educ6_5 = lin1b$coefficients[20]
  )

  ## Dynamic estimator
  rezzo <- 
    nlm(fedoc_ll,
        c(numeric(K), 0, seq(from=0.1,to=1,length.out = J-2)),
        k = k,
        y0s,y1s,y2s,y3s,
        X1s,X2s,X3s,
        h = 1
    )

  (r_fedol <- tibble(
    estimator = "FEDOL",
    h = 1,
    s = s,
    point = ifelse(s==0,"point","bs_rep"),
    b_ltdhi = rezzo$estimate[1],
    b_child = rezzo$estimate[2],
    b_married = rezzo$estimate[3],
    b_emp_2 = rezzo$estimate[4],
    b_emp_3 = rezzo$estimate[5],
    b_emp_4 = rezzo$estimate[6],

    rho = rezzo$estimate[7],

    gamma3 = rezzo$estimate[8],
    gamma4 = rezzo$estimate[9]
  ))

  # Do for a few different values of h.
  for(hh in c(0.1,10)){
    rezzo <- 
      nlm(fedoc_ll,
          c(numeric(K), 0, seq(from=0.1,to=1,length.out = J-2)),
          k = k,
          y0s,y1s,y2s,y3s,
          X1s,X2s,X3s,
          h = hh
      )

    (r1_add <- tibble(
      estimator = "FEDOL",
      h = hh,
      s = s,
      point = ifelse(s==0,"point","bs_rep"),
      b_ltdhi = rezzo$estimate[1],
      b_child = rezzo$estimate[2],
      b_married = rezzo$estimate[3],
      b_emp_2 = rezzo$estimate[4],
      b_emp_3 = rezzo$estimate[5],
      b_emp_4 = rezzo$estimate[6],

      rho = rezzo$estimate[7],

      gamma3 = rezzo$estimate[8],
      gamma4 = rezzo$estimate[9]
    ))

    r_fedol <- bind_rows(r_fedol,r1_add)
  }

  ## Bind them all
  rows_out <- 
      bind_rows(
        r_fedol,
        r_lin1b
      )

  return(rows_out)
}
